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Abstract. In disordered, viscoelastic or viscoplastic materials a sample response 
exhibits a recovery phenomenon after the removal of a constant load or after creep. 
We study experimentally the recovery in paper, a quasi two-dimensional system with 
intrinsic structural disorder. The deformation is measured by using the digital image 
correlation (DIC) method. By the DIC we obtain accurate displacement data and 
the spatial fields of deformation and recovered strains. The averaged results are first 
compared to several heuristic models for in particular viscoelastic polymer materials. 
The most important experimental quantity is the permanent creep strain, and we 
analyze whether it is non-zero by fitting the empirical models of viscoelasticity. We 
then present in more detail the spatial recovery behavior results from DIC, and show 
that they indicate a power-law -type relaxation. We outline results on sample-to- 
sample variation and collective, spatial fluctuations in the recovery behaviour. An 
interpretation is provided of the relaxation in the general context of glassy, interacting 
systems with barriers. 
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1. Introduction 

The nonlinear behavior of materials is highly topical thanks to several observations of 
"universality", that is generic behavior that is not dependent on the particulars of the 
example at hand. Examples for such generality can be found from the crackling noise or 
intermittent response by which many systems - magnets (Barkhausen noise), dislocation 
assemblies (strain avalanches), foams, granular systems, fracture and creep deformation 
etc. - react to external forces or influences [1] EJ EJ [H El E] . Likewise, universal behavior 
is found in "aging" and "rejuvenation", two concepts first introduced in the theory of 
spin glasses where at the simplest these can be taken to imply that the correlation 
or response properties of the system are found to be dependent on the observation or 
waiting time scale [7]. This dependency is moreover often non-linear. 

The interest in these questions has many origins. On one hand, one can look 
for classes of behavior - as in crackling noise or in the aging scenarios, where in 
both cases simple models lead the way by giving examples of distinct behaviors. On 
other hand, for a materials scientist or an engineer it is important to gain access to 
reliable effective descriptions of materials or "laws" and equations that summarize the 
nonlinearities. In particular, for viscoelastic and viscoplastic polymers this has meant 
the accumulation of vast quantities of semi-empirical viscoelastic theories and associated 
measurement results [E]. For creep, or material deformation under a constant load, there 
are also similar empirical bodies of data and "master equations" to explain their generic 
features a classical example being logarithmic creep deformation. There have been 
recent advances in treating "interacting elastic systems" in such creep conditions in 
statistical mechanics, in the presence of disorder ([9] and see also [10] for a connection 
to viscoplastic deformation). 

To combine the fundamental physical issues and practical (but often much more 
extensive) observations is one challenge for the statistical mechanics of materials 
deformation. In this work, we study the recovery of creep strain in paper in this 
context. Fiber networks (as paper) consist of a "frozen" random structure which 
changes only as a result of microscopic fracture damage and energy dissipation related to 
irreversible plastic deformation [11]. On the smaller scale, fibers are viscoelastic-plastic 
" beams" that themselves exhibit creep behavior if tested individually. The fact that the 
typical creep behavior of the network does not equal that of single fibers highlights the 
importance of coarse-grained (or "collective") behavior [HI [12]. Studies of the physics of 
paper have noted that paper has the interesting rheological property of delayed recovery. 
Following a stress-strain cycle or creep load, the remaining stress at zero strain is time- 
dependent ([LT1[13]). Similar mechanical phenomena exist in other "fibrous networks" 
that have received recently attention such as bucky- and graphene oxide paper and 
actine networks [El [151 HH]- In polymeric systems it has been suggested that the 
rheology is controlled by the motion of defects [17] . Extremely long term creep-recovery 
studies have been also conducted for polymers [IB] . In the case of crystalline materials, 
polycrystalline ice has been demonstrated upon unloading to show logarithmic delayed 
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deformation [19]. This is fairly close to what we find below. 

Here, we study the dynamics of strain recovery after the loading through one creep 
cycle. A sample is in each test creep-strained at a constant stress and then the strain 
decrease and the deformation field is measured as a function of time. A typical creep- 
recovery curve is presented in Figure [U In the creep-recovery the sample experiences 
first an initial recovery and after that delayed recovery follows. In the delayed recovery 
the rate of recovery decreases until the sample has recovered all the possible deformation 
(in paper physics, for a discussion of these issues see [201 EI])- The essential questions 
are thus two: i) what is the non-recoverable strain if any, and ii) what is the dependence 
of the relaxation on time? A further question is what can be stated of microscopic 
dynamics and its relation to the sample-level, coarse-grained response. 

By the usage of the Digital Image Correlation -technique (DIC) we can extract with 
very high accuracy the global sample strain and also the spatial strain fields, thus also 
the details of the creep and recovery phases [221 E3]. The main results that we gain 
are the following (Figure [U presents the definitions of the relevant quantities). First 
of all, by careful measurements by almost five orders of magnitude in time, we gain 
information on the ensemble- averaged behavior of the recovery. This can be contrasted 
with empirical models of viscoelasticity, and we also find that part of the creep strain 
remains permanent. Certain models work fairly in comparison with the data, while most 
are clearly not valid. It appears that the relaxation can be best fitted by a power-law 
with a small exponent, with - as stated already - a remnanent strain in the limit of infinite 
times. We study the sample-to-sample and spatial variations of the initial creep part and 
the corresponding recovery dynamics. The spatial recovery rates follow a distribution, 
whose width scales with time in the same manner as the average recovery rate. This 
indicates an asymmetry compared to the dynamics under loading, which means that 
there is no simple "superposition principle" relating the two dynamics under loading 
and recovery. 

The structure of this paper is as follows. In Section 2, we present the experimental 
details as far as the procedures and the DIC technique are concerned. In Section 3, we 
show the experimental data, and discuss some heuristic models and their comparison 
with the data. We also present a possible interpretation of the data via an energy 
landscape picture. This allows also to consider the role of internal stresses in the 
relaxation. Section 4 finishes the paper with a discussion. 

2. Methods 

2.1. Experimental setup 

Pictures are taken of the sample during the creep and the recovery process and the 
displacement field is measured by using the digital image correlation method. The 
setup is shown in Figure EJ A 100 mm long and 30 mm wide sample is attached to the 
upper and lower clamps. The upper clamp is attached directly to the frame while the 
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Figure 1. A schematic strain-time curve of a creep-recovery experiment with the 
initial strain, the delayed strain, the initial recovery and the delayed recovery marked 
in. The sum of the initial strain and the delayed strain is called the total strain and 
the sum of the initial recovery and the delayed recovery is called the total recovered 
strain. 

lower clamp with a mass of 855 g is let hang freely. The load is attached to the lower 
clamp, which can be moved to up and down positions using a pneumatic cylinder. The 
speed of the lift movement is set to 1 cm/s, so that the load is applied steadily. In its 
unloaded state the paper has an 1 cm slack, and the load is fully applied within order 
of one second after the lift movement has started. 

A global displacement measurement device, Omron laser distance sensor, is 
attached to the measurement device to follow the movement of the lower clamp. The 
results of the global displacement measurement is used for the consistency checking of 
DIC results, since its accuracy was found to be 0.1 mm due to mechanical limitations 
and the accuracy of the sensor. 

The sample was imaged with PCO's 1 mega pixel grayscale digital camera, 
SensiCam 370 KL 0562. This camera has a very low thermal noise ratio. An example 
of the images taken is in Figure Images were taken at an exponentially decreasing 
rate during the recovery: 18 images during the first hour and 18 pictures after one hour 
of recovery. 

The lower clamp has a mass of 855 g, and since the paper thickness is 0.1mm, 
the stress applied to the sample during imaging was 2.8 MPa. The sample was loaded 
fo 30 to 60 seconds during the recovery phase due to the imaging. The total elastic 
strain due to loading was 0.13 ± 0.01% since an elastic modulus 2.2 ± 0.1 GPa was 
measured. Strain due to lower clamp mass was well below yield strain 0.3% and the 
order 0.1% which is regarded as a linear viscoelastic limit of paper in engineering: in 
Ref. [23] it is shown that if the stress applied to paper is below 0.1%, nearly all of the 
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Figure 2. The setup for creep experiments. A load is attached to the lower clamp 
and its movement is controlled by using pneumatic cylinders. During the unloading 
phase the pneumatic cylinders let the sample recover freely without any load. When 
pictures are taken during the recovery, the load due to the weight of the lower clamp, 
is temporarily applied to the sample. 



Figure 3. A fibre scale image from paper. The image size is approximately 4x3 mm 
and there are about 950x700 pixels. The image is taken from the centre of the sample. 



deformation caused by the stress is recovered immediately after the stress is removed. 
The relative strain caused by the load application during imaging is of the order of error 
bars presented in the Fig. [6] and the disturbances due to the load application to the 
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displacement measurement are smaller. 

2.2. Digital Image Correlation 

The digital image correlation is the task of finding a deformation function, mapping 
coordinates from a reference image to coordinates in the test image [251 1221 |26] . The 
deformation function is presented as cubic splines, where knots are defined in an 
evenly spaced grid, with knot spacing hxh pixels (crate). The exact algorithm for 
the deformation computation is described in [27J . A cubic spline approximation of the 
deformation function leads to a locally minimized elastic energy and it can represent 
global affine deformations correctly. The method lies between "global" and "local" 
methods. In the global approach one defines a single criterion, e.g. affine deformation, 
which is globally optimized. In the case of local methods, one minimizes the error in 
each zone of interest, and a global criterion for the deformation function is not defined. 

In our analysis a crate size of 64 pixels was used. The boundary errors were 
minimized by printing black frames to the sample, and aiming the camera to the centre 
of the frames. The black colour printed with a regular office printer was uniform enough 
to prevent the algorithm from seeing any texture in it. 

The crate 64x64 pixels was chosen based on tests of the algorithm. For example, 
an artificial deformation resembling the experimental one was made to the image. The 
displacement fields computed with crate sizes of 16 and 32 pixels are in Figures H] and 
El respectively. With both crate sizes, the linear displacement is clearly seen but in the 
Fig. |5l crate size of 32 pixels, the displacement field is much smoother than in the Figure 
H] where a crate size of 16 pixels was used. The artificial deformation was found and 
a correct amount of deformation was obtained, but in the area where the deformation 
should be linear, an average error of about 0.2 pixels is seen. After a further change of 
the crate size the noise created by the algorithm is reduced significantly, being less than 
0.1 pixels. Increasing the crate size improves the accuracy, but reduces the details seen 
in the spatial displacement field. In order to measure the strain accurately, we used a 
rather large 64x64 pixel crate size, since the main purpose was to measure the global 
strain. In Ref. [28] there is additional discussion about the accuracy of the method. 

2.3. Strain measurement 

We measure the strain in experiments using two different methods: the global 
displacement on the image area and the spatial average of spatial strain rates. When the 
global displacement of the image area is used, we take an initial picture as a reference 
image and compare all consecutive images to the initial one. The global displacement 
between initial and the reference is measured using equation 

d=(y'-y) x , (1) 

where global displacement on the image area d is averaged over the x positions. 
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Figure 4. Image deformed 
with artificially made lin- 
ear deformation and a small 
bump. The crate size of the 
computation is 16 pixels. 
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Figure 5. Image deformed 
with artificially made lin- 
ear deformation and a small 
bump. The crate size of the 
computation is 32 pixels. 



For the spatial average of the local strain rate we define an evenly spaced grid 
on an image, which consists of a discrete set of points {(xx,yi), (^2,2/2)) •••}• The test 
and reference images are defined as consecutive images taken during the experiment. A 
spatial strain rate in a grid point is computed from 

& - Ail ( 2 ) 

€ ~ AtV {Z) 
where Ay is defined as an y-directional displacement of grid point (i, j). At is the time 
difference between images. The distance to the lower clamp is computed from fitting 
the plane to displacements and estimating the position where the displacement is zero. 
The spatial average of the strain rate is computed over the grid points, to then yield 
the recovery strain rate during the experiment. 



2.4- Measurements 

We performed 32 first creep first recovery experiments. The creep times in the 
experiments were mainly between 20 and 500 seconds in length. In two experiments 
the creep time was longer: 1230 and 3078 seconds. 

The creep in paper, delayed deformation depicted in Fig. (TJ can be divided, as 
for many materials, into three stages in time, primary, secondary and tertiary creep 
[2~T] . When the load is applied to the paper the primary creep (delayed strain) starts 
immediately. In the primary creep the strain follows Andrade's power law at relatively 
high stresses 

e(t) =at b + c, (3) 

where a, b and c are constants and t is the creep time [20]. In paper, the rate of creep 
does not exhibit secondary regime where the creep rate is constant, but rather the 
creep rate continues to decrease until the tertiary regime [2T]. Rheological phases are 
sometimes defined so that the deformation gained in the primary creep is recoverable, 
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and the secondary creep deformation is nonrecoverable [2U1EI]- Tertiary creep begins 
roughly when the strain rate starts to increase after reaching its minimum at the end 
of the secondary creep. After that, the strain rate accelerates until the sample fails. 

We chose a wide range of total creep strains as a starting point to the recovery. 
We observed a primary creep close to Andrade's law varying from sample to sample as 
e(i) ~ go that the exponent b was smaller than the 6=1/3 associated typically 

to Andrade creep. A sharp transition to secondary creep was not observed, but the 
exponent values b had a tendency to increase slowly during the experiment. 

The sample material used in the experiments was ordinary copy paper with a basis 
weight of 80 g/m 2 . As usual, industrial paper has an anisotropy which is denoted by the 
so-called Machine Direction (MD) and Cross Direction (CD): if strained in the principal 
directions they differ e.g. in the typical degree of ductility in a tensile test, where the 
CD turns to be much more ductile. 

In the creep phase of every experiment a load with a mass of 4071 g was applied to 
the sample in the CD. With this load we get an initial stress of 13.3 MPa. Environmental 
conditions were kept at constant 33 ± 1 % of relative humidity and 33 ± 1 °C of 
temperature. With the parameters used, the absolute humidity of the environment 
is equal to the situation, where the relative humidity is 50 % and temperature 23 °C, 
which is the standard operating point of mechanical testing of paper. The studies of 
Ketoja et al. showed that the master curve analysis - for logarithmic creep - can be 
used by introducing a moisture -dependent shift coefficient. They observed for the creep 
rate an exponential dependence, similarly to the elastic modulus of paper [29]. The fact 
that the dependence of the viscoplastic deformation on environmental conditions can be 
"scaled away" makes us expect that qualitatively the relaxation would not change as a 
function of the moisture. 

3. Results 

3.1. Overview 

In Figure [6] there is a typical strain curve from one creep- recovery experiment. The 
strain was measured approximately 3 seconds after the load was applied. The load was 
applied to the sample at t = and the sample responded with a strain of approximately 
0.55 % which consists of instantaenous initial strain and creep which starts to accumulate 
immediately when the load is applied to the sample. 

The initial strain is followed by a delayed strain that increases until the load is 
removed. In the measurement shown in Figure [6] the load was removed after 150 seconds 
of creep. Immediately after the removal, the sample underwent initial recovery which 
was estimated by measuring the strain approximately 3 seconds after the removal of 
the load. In Figure [6] the initial recovery is roughly 0.7%. After the initial recovery 
the recovery process of the sample continues, and the part taking place after the initial 
recovery is called the delayed recovery. 
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Figure 6. A typical example of the deformation curve in a creep-recovery experiment. 
The total strain is computed as the total displacement of the image area to the strain 
direction using a displacement field obtained from DIC using Eq. (TTJ). Error bars 
indicate a 0.1 subpixel limit. Pixel size is 4 /im. In the creep regime the sample 
deforms according to Andrade's law with an exponent e(t) ~ t . 



3.2. Averaged creep recovery 

Next we will present the averaged time- dependent recovery, and compare it to a number 
of heuristic models. These models start with the typical linear viscoelasticity-like 
Maxwell-Kelvin one, which consists of simple elements describing an ideal viscous and 
an ideal elastic material. It has the form 

e(t) = ae~i + d, (4) 

where a, b and d can be seen as independent fitting parameters. The parameter d 
expresses the amount of unrecoverable strain. Similarly the sum of the parameters a 
and d is the amount of the remaining strain immediately after the initial recovery. 

A stretched exponential model (also motivated by the Fancey model of viscoelastic 
response, [311 [32]) is also used, 

e(t) = ae~(*y + d, (5) 

where a, b, c and d are independent fitting parameters. Similarly to the Maxwell-Kelvin 
model also here the last parameter d relates to the amount of unrecoverable strain and 
the sum of a and d should be the remaining strain immediately after the initial recovery. 

Another heuristic model is Schapery's thermodynamical model, where stresses were 
independent state variables and the entropy production and the Gibbs free energy were 
specially defined [33] : 
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a [%} 


b 


c 


d [%} 


US 


SAIR 


Maxwell-Kelvin 


0.29 


1830 s 




0.90 


0.90 


1.19 


Fancey 


1.46 


5.82 s 


0.10 ± 0.02 


0.79 


0.79 


2.25 


Schapery 


1.61 


3730 


0.96 ± 0.01 






1.61 


Power law 


2.19 


0.11 ± 0.02 




0.58 


0.58 





Table 1. Fitting parameters for different models. The amount of unrecoverable strain 
(US) and strain after initial recovery (SAIR) obtained from the fitting parameters are 
in the rightmost columns. 

where a, b and c are independent fitting parameters. The parameter t is the creep time 
of the sample which in this case was 592 seconds. The parameter a tells us again the 
amount of strain left immediately after the initial recovery. 

Finally, we also tried the combination of a power-law decay and a constant, 

e(t) = (at)- b + d, (7) 

where a, b and d are independent fitting parameters. The parameter d is again the 
unrecoverable strain. 

In Figure [7] the four cases are tried with fits to the experimental recovery data. The 
data is computed using global displacement on the image area using Eq. jl]). 

The first, the Maxwell-Kelvin model is clearly proven insufficient. Fancey 's 
stretched exponential model is slightly better than Schapery's thermodynamical model, 
and finally the power law fit is at least as good as Fancey's model. However, a decisive 
conclusion between an exponential or a power law is difficult. 

The fitting parameters for the cases in Figure are shown in Table [U From the 
Table we conclude that the unrecoverable strain does not exhibit as much variation from 
model to model as the strain after initial recovery. However the unrecoverable strain 
values from the various fits are still quite different. Therefore, much longer experiments 
would have been necessary to estimate the real asymptotic behavior of the recovery. To 
half the recoverable strain would take roughly 500 times longer, according to the Fancey- 
style fit (Eq. (jSJ)). The power-law fit exponent (0.11) could be taken to correspond to a 
logarithmic decay, as an analogy to logarithmic creep. However, the role of fluctuations 
becomes in logarithmic creep more important with time, in contrast to what happens 
in relaxation which we discuss below. 

Only three parameters in Table [1] have their errors estimates presented because 
only those parameters are remaining nearly constant for all of the individual recovery 
curves. All of the other parameters are more or less related to the sample properties 
which vary from sample to sample. In other words the three parameters relate to the 
general shape of the recovery curve. 
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Figure 7. Averaged recovery as a function of the recovery time. Four model 
expressions fitted to the experimental recovery data. The recovery is obtained by 
comparing the image taken right after the initial recovery to later images using Eq. 
(UJ. Experimental data is average over 32 experiments. 



3. 3. Spatial deformation fields 

The process of delayed recovery is illustrated in Fig. |8] The displacement fields in these 
figures are calculated by comparing the image taken right after the initial recovery to 
the later images. The displacement is presented in micrometers. The total amount 
of spatial, recovered strain can be read from the colorbar next to each displacement 
field. From Fig. [8] it can be seen that the displacement field is heterogeneous during 
the delayed recovery. Such local dynamics in the displacement field continue upto the 
end of the experiment. The heterogenity of the displacement field is obvious, but the 
experimental setup does not allow us to draw very detailed conclusions for two reasons. 
The first is that our image area is limited to only a small part of the recovering sample, 
and the second reason is that we are integrating recovery over timescales which differ 
over several magnitudes. Thus we limit the consideration of the spatial deformation 
field to the essential properties of strain rate distributions and compare those to an 
analogous study done for primary creep [30] . 

In Figure [9] the mean and the standard deviation of spatial deformation rates are 
presented. The local deformation rate is computed according to the Eq. (j2J). The 
standard deviation Ae(t) is defined by 



Ai(t) 



\ 



1 

N 



(8) 



h3 
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Figure 8. The displacement fields after 30, 1700 and 823000 seconds of delayed 
recovery, respectively. Total amount of a local recovered strain is shown in micrometers 
in the color-bar. 



The local deformation rate is evaluated at 100x73 grid points. The decay of the 
averaged local deformation rate is equal to the recovery measured from the global 
displacement on the image area. 

The deviation from a power-law behavior at large times is probably attributed 
to the inability to measure very small strain rates and a subsequent biasing of the 
average over experiments. However, a possible physical explanation for the deviation 
cannot be ruled out. The generic features of the figure support strongly the above- 
mentioned observation, that a power-law relaxation fits the data best and extends over 
three decades: this then in turn should imply that the viscoelastic models discussed 
above are not valid. 

In Figure [10] we present normalized distributions of local deformation rates. The 
deformation rate from each experiment is normalized according to 



where m is the average deformation rate and s is the standard deviation of the spatial 
deformations in the experiment. Clearly, both the average recovery rate and its 
variations follow the same scaling in time. In Figure [11] we present as an example 
normalized recovery rate distributions in one experiment. The form of the recovery rate 
distribution remains unchanged during the recovery of the sample after a normalization 
with s as in Eq. (jUJ). Note that contrary to the previous Figure [TO] the data is not 
presented in a semi-log scale. 

For the creep recovery, the standard deviation Ae of the spatial strain rates and 
the spatial average (e) obey the same scaling in time Ae ~ (e) ~ t~ hl . In another study 
by us, on spatial and global strain rates in creep, the fluctuations was shown to exhibit 
a relative increase during Andrade's and logarithmic creep [30]. The Andrade's law was 
observed for the primary creep as (e) ~ i~ 0,7 , but the spatial fluctuations decreased 
during the primary creep according to Ae ~ t -0 55 . An analogous behavior of the ratio 
of creep rate fluctuations to the average creep rate was found for logarithmic creep, 
which is observed at lower stresses. The recovery process does not behave similarly, 
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Figure 9. Averaged recovery rate and the standard deviation of spatial recovery rates 
as a function of recovery time. The data is averaged over 32 experiments. The recovery 
rate agrees with the global recovery presented in the Fig. [7] The standard deviation 
shows that the relative width of the recovery rate distribution does not change during 
the recovery process. 



as the ratio of fluctuations to the mean recovery rate is constant. Thus, despite the 
fact that we observe a power law recovery the results suggest strongly that recovery is 
not "inverse Andrade creep" nor inverse logarithmic creep at the microscopic level. A 
direct consquence is that the applications of superposition ideas, familiar in the theory of 
models of viscoelastic response (see e.g. [18] ) should not work as the role of fluctuations 
is different in creep and relaxation. 



3.4- Statistical properties 

In this section we show briefly the sample to sample variation of the recovery process. 
Partially elastic behavior of paper can be observed from Figure [12] where the initial 
recovery is shown to correlate with the initial strain. Here we neglect the dynamics 
of the creep process during the first three seconds during the initial strain and initial 
recovery. The dashed line in Figure [12] represents ideal elasticity which means that the 
initial recovery is equal to the initial strain. We see that the correlation between the 
initial recovery and the initial strain is rather linear after the initial strain is over 0.6%. 
By estimating where the measured data approximated by a straight line intersects the 
line of linear elasticity, one concludes that the ideal elastic behavior takes place when 
the initial strain is less than 0.5 %. One can also consider the recovered strain at the 
end of the experiment as a function of the maximum strain at the point the creep phase 
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normalized strain rate [arb. units] 



Figure 10. Normalized recovery rate distributions averaged over 32 experiments. 
Comparison of three different recovery times and distribution is normalized according 
to Eq. ©. We see that the average recovery rate and its variations follow the same 
scaling in time. 




normalized strain rate [arb. units] 



Figure 11. Normalized recovery rate distributions from a single experiment. The 
distribution is normalized according to Eq. We note that the normalized form of 
the recovery rate distribution is unchanged during the experiment. 



is finished. It shows similar features as Figure [12) a roughly linear behavior for total 
strains over 0.6 %, with a slope less than one. 
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Figure 12. Initial recovery as a function of initial strain. Dashed line represents 
behavior of an ideally elastic material. 



3.5. Fitting a landscape model 

Next we consider an energy landscape model to look for interpretations of the relaxation 
result. The main observations to be utilized are that we know that the strain rate scales 
close to dtt ~ —1/t, and that the local fluctuations do not seem to become in relative 
terms more important with time. This means that we may try a mean-field model as 
the role of fluctuations will not increase as time progresses in the relaxation process. 
Note that coarse-grained models of crystal plasticity [38] have difficulties in reproducing 
creep and relaxation phenomena, though they include explicit elastic interactions of the 
local yield strains. 

Let us start with the general note that the local strain can be written as follows: 
e(x, t) = € p i(x) + e e ^ comp (x) +e e ij ree (x, t). Here, we state that at the end the load history 
will leave the material with a local irreversible strain (field) e p i. In equilibrium (if ever 
reachable), this results in internal stresses that need to be balanced, and this creates 
corresponding "as such" elastic strains to maintain the compatibility of local internal 
stresses: e e i )Comp [31]. In addition, we have the leftover, the free elastic strain e e ij ree 
which is a function of space and time, and it is the one that actually relaxes away and is 
being measured in the recovery. Let us consider the dynamics of the free elastic strain 
in the averaged sense and drop the dependence thereof on x. This means that in the 
absence of a theory which accounts correctly for the spatial interactions of the stress or 
strain components, we resort to a mean-field description, where such complexity from 
local dynamics is summarized with an averaged "field". 

Now, the driving process in the relaxation is the fact that the free elastic strain 
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is not distributed equally, and in the optimal way: the elastic energy content averaged 

Ee 2 

over the sample (-^f L ) x is higher than in the equilibrium. So what we can write now 
is an equation, that states that "strain change rate = derivative of elastic energy with 
respect to strain times a typical relaxation rate". This is dte = —Ee/r. This of course 
would indicate a simple exponential decay. The usual trick applied in other similar 
systems with relaxation phenomena (as magnetic ones and liquid crystals (35]) is to add 
a typical barrier U(e), which governs the relaxation at a given time. Thus we assume 
that such typical barriers at e dictate the rate, and that they change with time. Then 
the equation reads 

d t e = -Ee/r exp (-U(e)/kT) (10) 

as it is usually written. kT has the usual meaning. 

In our case we know that the left-hand side scales as —1/t. Since e is the 
corresponding time derivative times t plus a constant, we can just extract U (times 
1/kT) by dividing by the strain and taking a logarithm. A constant shift factor remains 
(E/t) which we set to zero for simplicity. We may allow for generality that e ~ A+Bt' 13 , 
whence 

-^w^lf^ <»> 

The important thing here is that for a reasonable 0, close to zero, the RHS goes as 
In (C/t), where C is a constant. If we now interpret t as strain solving it from the 
relation of time and strain, we have in principle reconstructed the distribution of the 
typical time-dependent and rate-setting barriers in terms of strain, as given by the model 
fit. One can plot the distribution and observe that such barriers get higher as the strain 
relaxes, as is natural given the slowing-down of the experimentally observed relaxation, 
compared to a simple exponential decay. This is illustrated in Figure [131 where we 
have utilized in the reconstruction of U(e) the parameters from the power-law fit, Eq. 
(I7j) from Table I. The schematic model presented here does not take into account in 
any direct way the local frustration |36j and localization of relaxation dynamics - if 
an unrelaxed area is surrounded by relaxed material, this will induce barriers to local 
rearrangements. The potential U should however be interpreted as a measure thereof. 

The next step is to interpret the relaxation, expressed as Eq. ( flQj) as a balance 
equation involving the elastic stress a and internal stresses a int is 

rd t e = -(a - a mt ). (12) 

The main point of the ansatz is that the relaxation stops once the internal stresses equal 
the tendency to relax. Assuming that a = Ee, one recovers the internal stresses. In 
terms of the barriers, they can be written utilizing Eq. ([TO]) again as 

(Tint = o"(l - exp(-U(e e i)/kT)), (13) 
which is equivalent to 

a int = Ee(l + ^) (14) 
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Figure 13. The result of Eq. (fTl]) for the effective landscape potential U(e). The 
constants have been chosen according to experiments, e = 1.3[%], and A = 0.58[%]. 



for the internal stress, e itself is time-dependent, and a correction appears proportional 
to r and the logarithmic time-derivative of the strain. 

Using the power-law fit for the strain, and parameters from the Table [TJ we can 
estimate the contribution of the elastic strain and the barrier term to the internal stress 
a int . One needs also to explore various guesses for r, the internal relaxation timescale. 
To illustrate the result, we plot the correction from the barriers to (Tint, or the expression 
y = —rd t e/Ee in the Figure [HJ The plot shows that if the relaxation time r is of the 
order of week up to years the main contribution of the relaxation of the internal stresses 
is due to the part proportional to the elastic stress. In other words, for reasonable 
values of r on time-scales up from one minute, y is always much smaller than unity and 
decays in an apparent power law fashion. This decay is in turn directly related to the 
relaxation exponent /3 indicated by experimental fit to a power-law form, and in this 
interpretation to the typical barrier dynamics. 

The simplest implication this results in is that there are three independent 
ingredients in the relaxation process: the inherent time-scale r, the particular dynamics 
of internal stresses that would be tightly coupled to the (relaxing) elastic ones, and the 
slowing-down due to barriers as such. 

4. Conclusions 

We have studied, in paper creep, the recovery of the creep strain. It turns out to be 
such that the phenomenon is characterized by two features: the statistical variations 
in each experiment and from sample to sample, and by the extremely slow response. 
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time [s] 



Figure 14. The expression y = —Td t e/Ee is shown using 1 week and 1 year relaxation 
times (t). The figure depicts the effect of the barriers to the relaxation dynamics by 
showing the magnitude of the barrier term in equation 1141 The constants have been 
chosen according to parameters from the Table [TJ 

The recovery seems to follow a power-law dependence on time. The local strain rates 
follow a probability distribution, which rescaled is independent of time. This indicates 
the presence of collective phenomena in the recovery dynamics. 

The main technique utilized here, the DIC, is clearly an advancement since it 
allows routinely to achieve good accuracy in the sample deformation measurement, and 
simultaneously gives access to the spatial strains. Thereby, the statistical behavior of 
local dynamics of relaxation in a single experiment and averaged over many experiments 
is available. 

We have presented empirical observation on the cross-correlations of the creep phase 
deformations and those of the relaxation phase. It concerns the immediate (elastic) 
relaxation and its relation to the creep strain at the beginning of the relaxation. This 
has all been done in one specific case, of material and creep test parameters, and it seems 
pertinent to underline that further work would seem important here. For instance, one 
should do cyclic experiments to mimic rejuvenation and aging studies in glassy systems 

BET]. 

The main results are perhaps, from the viewpoint of glassy systems, those that 
regard the recovery as a function of time and with the behavior of the distribution of 
local variations. It appears worthwhile to study the phenomena in greater detail. The 
results suggest that despite large variation in the initial creep deformation we see that 
the generic power-law recovery behavior remains the same. One feature one should pay 
attention to is the fact that the recovery fluctuations scale similarly to the mean recovery 



Relaxation of creep strain in paper 



19 



rate, in contrast to the original creep process where they increase in relative terms [30J. 
There, it appears that this is due to an underlying phase transition, controlled by the 
applied stress. 

The landscape fitting model gives way to interpret the slow dynamics via the 
effective barriers that dictate the relaxation rate. These are not included in viscoelastic 
models. Such a description makes more clear the observation that the relaxation is 
indeed of power-law form - contrary to any viscoelastic models. The derived typical 
barrier increases as the relaxation process advances, which might relate to the increasing 
resistance of the surrounding, more relaxed regions to local rearrangements. One can 
push the fitting to such a model further by making an Ansatz about the role of internal 
stresses, which leads to some understanding of the role of various mechanisms in the 
relaxation. It would be interesting and important to develop models that take into 
explicitly account the spatial part of the dynamics, and could thus be compared with 
such data as the DIC technique gives access to. One possibility might be extending 
mesoscopic models of plasticity in terms of elastic manifolds [38J . 
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